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Abstract 

In this paper we provide an extensive analysis of the global dynamics of high-area- 
to-mass ratios geosynchronous (GEO) space debris, applying a recent technique 
developed by Cincotta and Simo (2000), Mean Exponential Growth factor of Nearby 
Orbits (MEGND), which provides an efficient tool to investigate both regular and 
chaotic components of the phase space. 

We compute a stability atlas, for a large set of near-geosynchronous space debris, 
by numerically computing the MEGNO indicator, to provide an accurate understand- 
ing of the location of stable and unstable orbits as well as the timescale of their 
exponential divergence in case of chaotic motion. The results improve the analysis 
presented in Breiter et al. (2005) notably by considering the particular case of high- 
area-to-mass ratios space debris. The results indicate that chaotic orbits regions can 
be highly relevant, especially for very high area-to-mass ratios. 

We then provide some numerical investigations and an analytical theory that 
lead to a detailed understanding of the resonance structures appearing in the phase 
space. These analyses bring to the fore a relevant class of secondary resonances on 
both sides of the well-known pendulum-like pattern of geostationary objects, leading 
to a complex dynamics. 
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1 Introduction 



Recent optical surveys in high-altitude orbits, performed by the European 
Space Agency Im telescope on Tenerife (Canary islands), have discovered a 
new unexpected population of 10 cm sized space debris in near geosynchronous 
orbits (GEO). These objects sometimes present highly eccentric orbits with 
eccentricities as high as 0.55 (Schildknecht et al., 2004 and 2005). Follow- 
ing the initial guess of iLiou and Weaverl (120041 ) who suggested that this new 
population may be constitued by GEO objects with high area-to-mass ra- 
tios, recent nume rical and analytical investigations were performed to support 
this assumption (lAnselmo and Pardinl l2005l: iLiou and Weaveii . 120051). In ad- 
dition, these authors and others, such as IChaol fl2006h " and later IValk et al.l 



( l2008al ). presented some detailed results concerning the short- and long-term 
evolution of high area-to-mass ratios geosynchronous space debris subjected to 
direct solar radiation pressure. More specifically, these latter authors mainly 
focused their attention on the long-term variation of both the eccentricity and 
the inclination vector. Moreover, some studies concerning the effects of the 
Earth's shadowing effects on the motion of such space debris were given in 



Valk and Lemaitrd (l2008bl ). However, nobody ever dealt with the question to 



know whether these orbits are really predictable or not on the time scales of 
their investigations. 



The objective of this paper is twofold. The first goal is the investigation of the 
long-term stability of high area-to-mass ratio space debris subjected to direct 
solar radiation pressure, by means of the Mean Exponential Growth factor of 
Nearby Orbits (MEGNO) criterion. Second, still considering high area-to-mass 
ratios, we bring to the fore a relevant class of additional secondary structures 
appearing in the phase space. 



The paper is organized as follows. In Section [21 we focus our attention to 
the specification of the underlying model and we give some details about the 
numerical aspects of the method. In Section [31 for the sake of completeness, 
we dwell upon the detailed definition of the MEGNO indicator, also providing 
a review of its main properties, in order to understand the behavior of the 
chaos indicator. Then in Section [H in the framework of the validation o f 
our implementation, we retrieve the results obtained by lBreiter et al.l (l2005l ). 
W e also discuss the si gnificance of the time of integration, recently reported 
by [Barrio et al.l (|2007[ ). In Section [3, we apply the MEGNO technique in order 
to give a insightful understanding of the stability of high area-to-mass ratio 
space debris. More specifically, we show that the orbits of such peculiar space 
debris are extremely sensitive to initial conditions, especially with respect to 
the mean longitude and the semi-major axis. Second, we perform extended 
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numerical analyses, showing that the related 2- dimensional phase space is 
dominated by chaotic regions, in particular when the area-to-mass ratio is 
large. In addition, we also provide some results presenting the importance of 
the initial eccentricity value in the appearance of chaotic regions. Finally, in 
Section [6], we present extensive numerical and analytical investigations of the 
additional patterns which will be identified as secondary resonances. 



2 The model 



For the purpose of our study, we consider the modehng of a space debris 
subjected to the influence of the Earth's gravity field, to both the gravitational 
perturbations of the Sun and the Moon as well as to the direct solar radiation 
pressure. As a consequence the differential system of equations governing the 
dynamics is given by 

r = flpot + 0,(1 + a© 4- a,p , 

where Opot is the acceleration induced by the Earth's gravity field, which can 
be expressed as the gradient of the following potential 

f/(r. A, 0) = V V I — ) V^{sm (t)){Cnm cos m\ + Snm sin mX) , (1) 



n=Qm=Q ^ 



where the quantities Cnm and Snm are the spherical harmonics coefficients 
of the geopotential. The Earth's gravity field adopted is the EGM96 model 



(ILemoine et al.l . 119871 ). In Eq. ([I]), /i is the gravitational constant of the Earth, 
Re is the Earth's equatorial radius and the quantities (r. A, 0) are the geocen- 
tric spherical coordinates of the space debris. are the well-known Legen- 
dre functions. It is worth noting that the potential of Eq. ([1]) is subsequently 
expressed in Cartesi an coordinates by means of the Cunningham algorithm 



( Cunningham . IQTOl ). 



Both the accelerations aq^ and a© result from the gravity interaction with a 
third body of mass m*, where * = 1 and * = 0, and can be expressed with 
respect to the Earth's center of mass as 




where r and are the geocentric coordinates of the space debris and of the 
mass m^, respectively. The quantity is the gravitational constant of the 
third-body. In our implementation, we chose the high accurate solar system 
ephemeris given by the Jet Propulsion Laboratory (JPL) to provide the posi- 
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tions of both the Sun and the Moon (IStandishl . Il998l ). 



Regarding direct solar radiation pressure, we assume an hypothetically spheri- 
cal space debris. The albedo of the Earth is ignored and the Earth's shadowing 
effects are not taken into account either. The acceleration induced by direct 
solar radiation pressure is given by 



A 



r© 



m ||r — TqI 



where Cr is the adimensional reflectivity coefficient (fixed to 1 further on in 
this paper) which depends on the optical properties of the space debris surface; 
Pr = 4.56 X 10^^ N/m^ is the radiation pressure for an object located at the 
distance of 1 AU; aQ = 1 AU is a constant parameter equal to the mean 
distance between the Sun and the Earth and Vq is the geocentric position 
of the Sun. Finally, the coefficient A/m is the so-called area-to-mass ratio 
where A and m are the effective cross-section and mass of the space debris, 
respectively. 



3 The Mean Exponential Growth factor of Nearby Orbits 



For the sake of clarity we present in this section the definition and some prop- 
erties of the MEGNO criterion. 



Let 7^(p, q), with p G M", q G T", be a n-degree of freedom Hamiltonian 
system and let us introduce the compact notation x = {p, q) G M^" as well 
as / = {—dTC/dq,dT-C/dp) G M^", then the dynamical system is described by 
the following set of ordinary differential equations 



d 
di 



x{t) = fix{t),cx) 



X G 



n2n 



(2) 



where o: is a vector of parameters entirely defined by the model. Let = 
(j)(t] XQ,to) be a solution of the flow defined in Eq. ([2]) with initial conditions 
(to, Xq), then it has associated the Lyapu nov Characteristic Number (hereafter 



LCN), defined by (IBenettin et all . Il98nh 



A = lim - In 



|<5</,(to 



(3) 



where S^{t), the so-called tangent vector, measures the evolution of an ini- 
tial infinitesimal deviation ^^(to) = between (j){t) and a nearby orbit, and 
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whose evolution is given by the variational equations (terms of order 0{S'^) 
are omitted) 

S^ = j^S^{t) = J{<p{t))d^{t), with J(0(t)) = |^(0(t)), (4) 

where J{(j){t)) is the Jacobian matrix of the differential system of equations, 
evaluated on the solution Let us note that the definition of LCN, given 
by Eq. ([3]), can also be written in an integral form 

A=limi r^-Mds, 



t^oo t Jo 5^{s) 

where Sa, = \\SJ\, 6^ = d^h- S 



The Mean Exponential Growth factor of Nearby Orbits Y^jt) is based on a 



modi fied time- weighted version of the integral form of LCN (jCincotta and Simd . 



20001 ) ■ More precisely 



t Jo d^{sj 

as well as its corresponding mean value, to get rid of the quasi-periodic oscil- 
lation possibly existing in Y^{t) 

1 /•* 



Y^{t) = - I Y^{s) ds 

b 



In the following we will omit the explicit dependence of Y and Y on the specific 
orbit 0, when this will be clear from the context. 

Actually, Y{t) allows to study the dynamics for long time scale s, where generi 



cally y (if:) does not converge, while lim^^oo Y{t) is well defined (ICincotta et al. 



20031 ). Consequently, the time evolution of Y{t) allows to derive the possible 
divergence of the norm of the tangent vector 8{t), giving a clear indication of 
the character of the different orbits. Indeed, for quasi-periodic (regular) or- 
bits, Y{t) oscillates around the value 2 with a linear growth of the separation 
between nearby orbits. On the other hand, for chaotic (irregular) motion, the 
norm of 6 grows e xponentially wi t h tini e, and Y{t) oscillates around a lin- 



ear divergence line. ICincotta et al.l (120031 ) showed that, for the quasi-periodic 



orbits, Y{t) always converges to 2, that is a fixed constant. Moreover, it has 
been shown that ordered motions with harmonic oscillations, i.e. orbits very 
close to a stable periodic orbit, tend asymptotically to Y(t) = 0. 



These latter properties can also be used to compute efficiently a good estima- 
tion of LCN, or similarly the Lyapunov time Tx = 1/A, by means of a linear 
least square fit of Y{t). Indeed, in the case of an irregular orbit, the time 
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evolution of Y{t) may be easily written as 



Y{t) ^ tti^t + d, t — i> oo, 



where a^, is simply related to LCN by the relation = A/2 and d is small for 
irregular and chaotic motion. Bur for regular orbits, after a transitory time, d 
is not necessarily close to zero. Thus, the value of d may be considered as the 
measure of how long the o rbits sticks to a regular torus before getting chaotic 
(ICincotta and Simol . I2OOOI ). 



Regarding the num erical computation of the MEGNO indicator, we adopt the 
same strategy as in iGozdziewski et al.l (I2OOII ). To be specific, in addition to 
the numerical integrations of both the equations of motion and the first order 
variation equations, we consider the two additional differential equations 



d 



6S 



d 



w 



(5) 



dt d ■ S dt 
which allow to derive the MEGNO indicators as 

Y{t)=2y{t)/t, Y{t)=w{t)/t. 

The MEGNO criterion, unlike the common Lyapunov variational methods, takes 
advantage of the whole dynamical information for the orbits and the evolution 
of its tangent vector, which results in shorter times of integration to achieve 
comp arable results. Moreover^ a couple of applications fo u nd in the literature 
'e^ 



Gozdziewski et al. 2001. 2008: Breiter et al. 2005: Cincotta and Simo 



2OOOI ) justify and confirm that MEGNO is relevant, reliable and provides an 
efficient way for the investigation of the dynamics by detecting regular as well 
as stochastic regimes. 



3.1 MEGNO and numerical integrations 



As previously mentioned, in order to evaluate the MEGNO indicator, we have 
to integrate the differential system of Eq. ([2]), the linear first order variational 
system of Eq. (jl]), as well as the two additional differential Eq. ([5]). We choose 
to write both the expressions of the perturbing forces and the variational sys- 
tem, i.e. the Jacobian matrix, in rectangular coordinates positions-velocities. 
In such a way we can overcome both the null eccentricity and the null inclina- 



tion singularity present in the dynamics of space debris (jValk et al.l . l2008d ) . 
Moreover, the explicit analytical expressions of the vector fields allow us to 
avoid the difficulties inherent in the classical method of neighboring trajecto- 
ries (two particles method). 

In order to integrate numerically the two differential systems of equations. 
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we adopted the variable step size Bulirsh-Stoer algorithm (see e.g. Bulirsh 
and Stoer, 1966, and Stoer and Bulirsch, 1980). Let us note that, for the 
purpose of validation, the numerical integrations were also made with a couple 
of other numerical integrators. However, the Bulirsh-Stoer algorithm seems to 
be the best compro mise between accuracy and efficiency. More over, as quoted 



bv IWisdomI (119831 ): What is more important for this study, Wenettin et al. 
1 198d) found that the maximum LCM^ did not depend on the precision of 



their calculation. It appears likely that as long as a certain minimum precision 
is kept, maximum LCE's may be accurately computed, even though it is not 
possible to precisely follow a specified trajectory for the required length of time. 

Although this latter observation was formulated in the framework of both Lya- 
punov variational method and Hamiltonian systems, it seems that it remains 
relevant in the computation of the MEGNO criterion, at least in the particular 
case of our analysis. 



3.2 Influence of the initial tangent vector do 



By constr uction MEGNO depends on the initial value of the tangent vector Sq as 
the LCE ( iBenettin et al.l . Il980l ). That's way we preferred to adopt the strat- 
egy of initialize randomly the initial tangent vectors in order to avoid some 
parts of the artificially created zones of low M EGNO due to the pro ximity of Sq 
to the minimum Lyapunov exponent directi on ( iBreiter et al.l . l2005l ) . Moreover, 
as pointed out by iGozdziewski et al.l (120011 ) , the random sampling of So is rel- 
evant in the sense that different initial tangent vectors can lead to different 
behaviors of the MEGNO time evolution while considering the same orbit. This 
observation has been reported in the framework of extra-solar planetary sys- 
tems and seems to remain similar in the case of Earth orbiting objects and 
more generally for high- dimensional dynamical systems (having more than 3 
degrees of freedom). 



Regarding the impact of the choice of the initial tangent vector 6o, we per- 
formed a set of exhaustive numerical investigations of regular orbits. More 
specifically, we compared the time-evolution of MEGNO using different initial 
tangent vectors and identical generic initial conditions. The results confirm 
that the random choice of the initial tangent vector induces a significant ran- 
dom behavior in the way MEGNO approaches the limit value 2, hence preventing 
this information from being useful to check the stability/instability character 
of regular orbits. Actually, when considering a slightly perturbed two-body 
problem (such as the central attraction disturbed by the oblateness of the 
Earth), the way MEGNO converges to 2 is completely unpredictable , leading to 



Lyapunov Characteristic Exponent. 
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more or less 50% of convergence of Y{t) to 2 from above and the other remain- 
ing 50% from below. This result is formally discussed in the following subsec- 
tions. However, when the order of magnitude of the perturbation is larger, the 
result does not completely hold anymore. In particular, when considering the 
perturbing effects induced by the 1:1 resonance, the MEGNO evolution no longer 
depends on the random choice of the initial tangent vector. In this latter case, 
the intrinsic stability of the chosen orbits se ems also to dictate the evolution of 
MEGNO as reported in ICincotta et al.l (120031 ). More specifically, the stability of 
the orbit seems to influence the time evolution of MEGNO the more the orbit is 
closer to a stable or unstable equilibrium point. For instance, regarding the or- 
bits extremely close to a stable equilibrium point, MEGNO generally approaches 
slowly the limit value 2 from below, even though some infrequent orbits present 
a MEGNO convergence from above. Conversely, the orbits initially close to the 
separatrices generally present a MEGNO approaching the value 2 from above. 



3.3 MEGNO for integrable systems 



In this section we will study the MEGNO indicator for integrable Hamiltonian 
systems and we will show that generically (if the system is not isochronous) 
it always converges to 2, moreover the way Y{t) reaches this limit value, say 
from higher or lower values, depends only on the choice of the initial tangent 
vector and not on the orbit itself. 



So let us consider an integrable Hamiltonian system and suppose to write 
it in action-angle variables, 7i = 7i(p), where p G -B C denotes the 
action variables and g G T" denotes the angle variables. Then the Hamiltonian 
equations are 



p = 0, 




The tangent space (to a given orbit) can be split into the action and angle 
direction, namely 6 = {dp, dg), thus the variational system can be written as 



Sp = 0, 

4 = = M{p) 5p . 



If the system is isochronous then M = 0, thus 6p and 6q are constant and 
Y{t) = for all t. On the other hand, if the system is non-isochronous we get 
5p{t) = <5p(0) and 5q{t) = (5^(0) + M(p(0)) 6p{0)t. To simplify the notations. 
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let us introduce 

M (p(0)) = Mo, (5p(0) = ^0 and 5,(0) = ryo • 
Using the definition of MEGNO, we get 



^'W - 7 / 



and tliis integral can be explicitly computed, obtaining 



s ds . 



Y{t) 



2 - ; log[l + 2Molo ■ + iMo^oYt'] + 



t (Molo)' 



arctan 



(Molo)^ - (Molo ■ r7o)^ 



One can check that the square root is well defined, i.e. positive, and thus one 
can cast Eq. into 

If L 

where Fi and F2 are positive functions and F2 is bounded. We can then con- 
clude that (see Fig. [T]) 



(1) if Mo^o ' Vq> ^ then Y{t) approaches 2 from below; 

(2) if Mq^o ■ '^0 < then Y{t) approaches 2 from above, in fact for large t 
the first contribution dominates the bounded term F2. 



In this last part we will consider if and under which assumptions the previous 
results concerning the convergence F — >■ 2 are still valid, for a quasi-integrable 
Hamiltonian system of the form H{p, q, e) = Hq{p) +eV{p, q). The main idea 
is the following: fix e > 0, but small, and consider a "non-chaotic" orbit (p^, 
namely an orbit without a positive Lyapunov exponent (or "if you prefer" with 
a bounded MEGNO) , then if e is sufficiently small this orbit is a perturbation of 
an orbit existing also for e = 0, 0o, and we can check that Y^^ = Y^^ + C(e), 
hence the smallness of such e-correction cannot change "the way Y goes to 2" . 
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T 




Time Time 



Fig. 1. MEGNO for quasi-integrable adimensional Hamiltonian system. We consider 
the evolution of 1^^ for the system H = p1/2+p2+ecosqi+ecos{qi — q2)- On the left 
panel e = 10~^, while on the right panel e = 10~^. In both cases e is small enough 
to confirm the theoretical predictions; let observe that in this case the matrix M is 
given by ( Q Q ) and thus the sign condition reads MSp^ ■ Sq^ = (5p 0*^9,0 • "^^^ °f 
time corresponding to 1/10 of period of the orbit. 



More precisely, the Hamilton equations are now 



. _ _dn _ _ dV 
^ dq dq 

^_dn _ ^ ^ dV 
dp dp ' 

and a similar decomposition can be done for the variational system 



^ ^ dpdq ^ ^ dq^ 

\ dp"^ dp"^ j ^ dpdq ' 

Looking for 5p and 8q as e-power series, i.e. 5p = 6pfi + edp^i + . . . and 
= ^q,o + + • • • , and collecting together, in the definition of MEGNO, 
terms contributing to the same power of e, we can thus get 



(Modp, 



Mo6pfi ■ Sg^o 



t Jo {Sp^o^ + {Sp,oy + 2Mo5p,o ■ dq^os + {Mo5qflfs^ 



sds + 0{t) 
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4 Validation of the method 



To validate our method we first apply the technique on a simplified model, 
containing only the Earth's gravity field expanded up to the second degree 
and order harmonics, namely, J2 = —C20, C22 and 5*22- For the purpose of the 
analysis, we followed a set of 12 600 orbits, propagated over a 30 years time 
span, that is th e order of 10^ fundamenta l periods (1 day) empirically r e quire d 
by the method (iGozdziewski et al.l . l200ll ). As reported in lBreiter et al.l (120051 ) . 
a 30 years time span seems to be relatively small for long-term investigations 
of geosynchronous space debris. However, the numerical integration of vari- 
ational equations in addition to the extrapolation of the orbit is quite time 
consuming. Indeed, the simulation with an entry level step size of 400 seconds 
takes approximately 20 seconds per orbit when including only the Earth's 
gravity field, whereas it takes 42 seconds with a complete model. Thus, the 
examination of large sets of initial conditions can take a lot of time (typically 
5 days for 10^ orbits). On the other hand, the analysis of the following section 
will bring to the fore some indications about the Lyapunov times (smaller 
than 30 years). As a consequence, the integration time can be considered as 
sufficiently large in the particular case of our study. 



For the purpose of this validation study, we consider a set of initial conditions 
defined by a mean longitude A grid of 1°, spanning 90° on both sides of the 
first stable equilibrium point and a semi-major axis a grid of 1 km, spanning 
the 42164 ± 35 km range. The other fixed initial conditions are eo = 0.002 for 
the eccentricity, zq = 0.004 rad for the inclination, = cuo = rad for the lon- 
gitude of the ascending node and the argument of perigee, respectively. These 
values have been fixed to compare o u r resu lts for the nearly-ge osynchronous 



orbits with the ones of iBreiter et al.l (120051 ) . As pointed out by iBreiter et al 



(I2OO5I ). due to the 1:1 resonance, good variables to present our results will be 
(oo, o"o), where oq is the osculating initial semi-major axis and a is the so-called 
resonant angle, i.e. a = \ — 9, where 9 is the sidereal time. 



Figure [2] (left panel) shows the MEGNO values computed using 30 years of inte- 
gration time. We identify clearly a blow-up of the typical double pendulum-like 
pattern related to the 1:1 resonance. Here, we plot only over a horizontal range 
of 180°, i.e. only one eye. The existence of both the stable and the two un- 
stable equilibrium points can be easily inferred. We observe that the phase 
space seems to be essentially filled in with MEGNO values Y(t) ^ 2, that is 
plenty of regular orbits. Moreover, the two separatrices are also identifiable 
and are associated with neighboring MEGNO values 2 < Y{t) < 4. Therefore, 
following the properties defined in Section [3], one could consider that these 
orbits are chaotic. However, we will show that this conclusion is false. Indeed, 
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Resonant angle [degree] Resonant angle [degree] 



Fig. 2. Tlie MEGNO computed as a function of initial mean longitude Aq and osculating 
semi-major axis qq. The equations of motion include the central body attraction as 
well as the second degree and order harmonics J2 , C22 and ^22 • The mean longitude 
grid is 1° and the semi-major axis grid is 1 km, spanning the 42164 it 35 km range. 
The initial conditions are cq = 0.002, iq = 0.004 rad, ^Iq = loq = rad. Time at 
epoch is 25 January 1991. The patterns have been obtained using two different times 
of integration, tf = 30 years [left] and tf = 300 years [right]. 

a careful identification of the MEGNO time evolution shows that the latter al- 
ways approach slowly the limit 2 from above. The closer to the separatrice, 
the slower the convergence. More precisely, orbits close to the separatrix in- 
tegrated over long time span present a bounded MEGNO evolution. Hence they 
should be considered as non chaotic. 



To clarify this point, we performed a similar study, but using a significantly 
longer time-span, namely 300 years. The results are shown in Figure [2] (right 
panel). For the sake of comparison, the color bars have been taken identical on 
both plots. Let us observe that the maximum value reached by the MEGNO is 4 
in the left panel and 2.5 in the right one. In the 300 years simulation (Figure [21 
right), the MEGNO values, associated with orbits close to the separatrices, turn 
out to be, on average, smaller than in Figure [2] (left panel), reaching almost the 
limit Y{t) — i> 2, due to the longer time of integration. Similarly, the dark zone 
in the neighborhood of the stable equilibrium point, corresponding to MEGNO 
values close to zero, is strongly shrunk, supporting the result that, in the limit 
of infinitely large t, only the orbit originating from the exact stable equilibrium 
point leads to Y = 0, whereas the neighboring trajectories converge slowly to 
F(t) = 2. 



Let us note that the iraportan ce of the integration time has been recently 
reported by lBarrio et al.l (120071 ) in the framework of applications of the MEGNO 
method. And, we confirm that a too short time of integration can give wrong 
conclusions about the dynamical behavior. Moreover, the latter paper also 
underlines some spurious structures appearing in the maps of the variational 
chaos indicators, explaining the presence of the sine wave of lower MEGNO with 
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a bulge at the center of Figures [21 ^''suggesting that the same periodic orbit 
is more or less regular depending on the initial conditions choice" . Actually, 
accordly to the latters authors and to our analysis, this conclusion is wrong 
because this spurious structure is related to numerical artifacts. 



5 High area-to-mass ratios analysis 



The study of the long-term stability of near-geosynchronous objects has re- 
cently prompted an increasing interest of the scientific community. In the 
particular case of classical near-geosynchronous objects, the problem has been 
solved by computing the MEGNO indicator for a family of simulated geosta- 
tionary, geosynchronous and super-geosynchronous orbits. A classical near- 
geosynchronous object has a period close to one sidereal day and is subjected 
to the main gravitational effects of the Earth, including the 1:1 resonance, 
luni-solar perturbing effects, as well as solar radiation pre ssure associat e d to a 
sma ll area-to-mass ratio (A /m <C 1 m^/kg). According to lBreiter et al.l (120051 ) 
andl Wytrzyszczak et al. ( 2007 ). the near-geostationary region presents chaotic 
orbits only very close to the separatrices, due to the irregular transits between 
the libration and the circulation regimes. Regarding the super-geostationary 
orbits, all of them seem to be entirely regular on the time scale of the inves- 
tigations, that is a few decades. 



The aim of this section is to provide a more extensive analysis of the dy- 
namics of near-geosynchronous space debris with high area-to-mass ratios 
{A/m ^ 1 m^/kg), subjected to direct solar radiation pressure. Our main 
objective is to study the effects of high area-to-mass ratios on the stability of 
the principal periodic orbits and on the chaotic components. This analysis is 
divided into three parts. First, in section 15.11 we focus our attention on the 
sensitivity to initial conditions; then, in section 15.21 we report the results of 
dedicated numerical analyses which emphasize the importance of the area-to- 
mass ratio value. Finally, in section 15.31 we study the influence of both the 
initial eccentricity and time at epoch. 

Let us recall that for large area-to-mass ratios [A/m > 10 m^/kg), the solar 
radiation pressure may become th e major perturbati on, by far larger than the 
dominant zonal gravity term J2 ( Valk et al. . 2008al ). In this particular case. 



the larger the area-to-mass ratio, the more affected the dynamics of the near- 
geosynchronous space debris, leading to daily high-amplitude oscillations of 
the semi-major axis, yearly oscillations of the eccentricity as well as long-term 
variations of the inclination. As an illustration. Figure [3] shows the orbital 
elements histories of the first 210 years of a geosynchronous high area-to-mass 
ratio space debris {A/m = 10 m^/kg). The yearly variation of the eccentricity 
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Fig. 3. Time-evolution of high area-to-mass ratio space debris. Orbital elements 
over 210 years for a A/m = 10 m^/kg; initial conditions are: oq = 42166.473 km, 
eo = 0.002, iQ = 0.004 rad, = ujq = rad and Mq = 4.928 rad. Time at epoch is 
25 January 1991. 
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reaches 0.2, which confirms the expected values predicted (e.g. Anselmo and 
Pardini, 2005, and Liou and Weaver, 2005). The inchnation evolution presents 
a well known long-term variation whose period is directly related to the area- 
to-mass ratio value. Regarding the longitude of ascending node as well as the 
argument of perigee, they both present a librat ion due to the cho sen set of 
initial condi tions. For further details, we refer to IValk et al.l (l2008al ) as well as 



Chad (120061 ) . where a full description of the long-term motion of high area-to- 



mass ratios space debris is given. 



5.1 Sensitivity to initial conditions 



To start with, we follow the evolution of two high area-to-mass ratio space 
debris {A/m = 10 m^/kg) defined by two sets of very close initial conditions, 
differing only in the 10th digits in mean longitude. Figure [3] shows the first 
one and Figure H] shows the second near orbit. We observe that the most dif- 
ference (in the behavior) take place in the semi-major axis and resonant angle 
panels. We notice that there are some differences in the dynamics of the semi- 
major axis already after 20 years and at the end of the integration. This is 
the same for resonant angle, confirming the hypothesis that the sensitivity to 
initial conditions is especially relevant for the semi-major axis and resonant 
angle whereas the difference between the other orbital elements remains small. 
We first focus our attention on the time evolution of the semi-major axis and 
resonant angle. As a complement to Figure [3l we numerically computed two 
orbits for two space debris with different area-to-mass ratios, A/m = 1 m^/kg 
and A/m = 10 m^/kg, whose initial conditions have been chosen near the 
separatrices, to emphasize their chaotic behaviors. 



Figure shows a blow-up of the evolution of the semi-major axis (top panels) 
and resonant angle (middle panels) over the time span of 250 years. It is 
clear that the semi-major axis presents some irregular components over its 
evolution, related to some transitions between different regimes of motion, 
clearly identifiable in the resonant angle plots. In addition, we also computed 
the corresponding MEGNO time evolution. The bottom panel in each graph 
shows the time evolution of the MEGNO indicator as well as its corresponding 
mean value. First, we see that the time evolution of Y{t) presents a quasi- 
linear growth almost since the beginning of the integration process, leading to 
the conclusion that these orbits are clearly chaotic over that time scale. 
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Fig. 4. Time-evolution of liigh area-to-mass ratio space debris. Orbital elements over 
210 years for a A/m = 10 m^/kg; the initial conditions are the same of Figure [3] 
but differ in the 10th digits in mean longitude. 

Therefore, we also computed the linear fit F (t) ~ t + c? in order to evaluate 
the Lyapunov time Tx. Tx is the inverse of the LCN (A) calculated by the 
linear regression coefficients = A/2. Let us remark that to avoid the initial 
transient state, the least square fits were performed on the last 85% of the 



16 



time interval. This latter analysis brings to the fore the fact that larger area-to- 
mass ratios lead to smaller Lyapunov times, i.e. larger Lyapunov characteristic 
numbers. Indeed, for A/m = 1 m^/kg, the Lyapunov time turns out to be 
on the order of 11 years, whereas it reaches the value Tx ~ 3.7 years for 
A/m = 10 m^/kg. 
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Fig. 5. For each graph, we show the orbital evolution of the semi-major axes a (solid 
line) superimposed with the evolution of the inclinations (dashed line) [top panels]. 
The time evolution of the resonant angles [middle panels] and the time evolution 
of the MEGNO indicator (Y and Y =< Y{t) >), as well as the corresponding linear 
fit Y{t) a^^t + d [bottom panels]. The area-to-mass ratios are A/m = 1 m^/kg in 
the upper panel and A/m = 10 m^/kg in the lower one. The initial conditions are 
chosen near the separatrices. The computed linear regression coefficients are given 
by = 0.043 (for A/m = 1 m^/kg) and = 0.134 (for A/m = 10 mVkg). 
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Second, let us also remark that the behavior of the MEGNO indicator is of 
particular interest in these cases. A careful analysis of Y{t) underlines some 
irregular patterns directly related to the evolution of a, in particular when the 
orbits seem to transit across the separatrices. Finally, we can also highlight the 
fact that the sudden changes between libration and circulation regimes occur 
mainly when the inclination changes its sign of variation, especially at the max- 
imum value for A/m » 1 m^/kg and at the minimum for A/m < 1 m^/kg 
(Figure O top panels, dashed line), with an empirical long-term periodicity of 
Tq, that is the long-term periodicity of the longi tude of the ascending node. 



which is all the more smaller when A/m is large (IValk et al.l . l2008al ) 



5.2 Extended numerical analyses 



We considered a set of 12 600 simulated orbits with various initial semi-major 
axes and mean longitudes. We took into account the following perturbing 
effects: second degree and order harmonics ( J2, C22 and S22), the luni-solar in- 
teraction as well as the perturbing effects of the solar radiation pressure with 
four values of the area-to-mass ratio {A/m = 1,5, 10,20 m^/kg). The results 
are reported in Figure O 



In the case with A/m = 1 m^/kg (top left panel) we recognize the same 
pendulum-like pattern as in Figure [2l Considering the same integration time 
(30 years), we notice that the MEGNO values tend to be slightly larger than 
in Figure [2] (left). Moreover, some irregularly distributed MEGNO values are 
clearly visible close to the two saddle uns table stationary poin ts. These results 



completely agree with those presented by iBreiter et al.l (120051 ) , where the solar 



radiation pressure was taken into account, but only for very small area-to- 
mass ratios (typically 0.005 m^/kg). Indeed, our latter analysis shows that in 
addition to the luni-solar perturbations, solar radiation pressure with small 
to moderate area-to-mass ratios, that is0<y4/m<l m^/kg, do not change 
considerably the phase space pattern. 

On the other hand, the remaining panels of Figure [H] show that the phase 
portrait becomes significantly more intricate with increasing area-to-mass ra- 
tios. Indeed, the width of the stochastic zone in the neighbourhood of the 
separatrices becomes relevant, with a large displacement of the separatrices 
on the phase plane. The larger chaotic region can readily be explained by the 
osculating motion of the separatrices due to the before-mentioned daily varia- 
tions of the semi-major axis with respect to some mean value as well as by the 
increasing amplitudes of the eccentricities. These variations lead inevitably to 
transits through both the regions separating libration and circulation motion 
for orbits initially close to the separatrices. 
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Fig. 6. The MEGND computed as a function of initial mean longitude Aq and ini- 
tial (osculating) semi-major axis ap. The equations of motion include the central 
body attraction, the second degree and order harmonics J2,C22 and 522, the lu- 
ni-solar interaction as well as the perturbing effects of the solar radiation pressure. 
The mean longitude grid is 1° and the semi-major axis grid is 1 km, spanning 
the 42164 it 35 km range. The initial conditions are cq = 0.002, iq = 0.004 rad 
and r^o = <^o = rad. The integration time is 30 years from epoch fixed at 25 Jan- 
uary 1991. The patterns have been obtained using four different area-to-mass ratios, 
A/m = 1, 5, 10, 20 m^/kg, represented, respectively in the top left, top right, bottom 
left and bottom right panel. 

Moreover, it is also clear that the usual double pendulum-like phase space 
shows a tendency to be distorted with an apparent displacement of the unsta- 
ble equilibrium points, whereas the stable equilibrium points remain almost 
fixed. This last result is however quite awkward insofar as there is no physi- 
cal interpretation to this phenomenon. Indeed, direct solar radiation pressure 
does not depend explicitly on the (mean) resonant angle with respect to the 
long-term investigations after averaging over short perioding terms. Therefore, 
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Fig. 7. Cartoon to illustrate the difference between mean and osculating initial con- 
ditions with respect to the semi-major axis (s.m.a.) evolution. For the sake of sim- 
plicity, the mean semi-major axis does not present any long-term variation, whereas 
the osculating semi-major axis presents daily oscillations related to direct solar ra- 
diation pressure (the implicit underlying model is radiation pressure only). It is 
clear that even if the osculating initial conditions af^^ and a2^^ are identical, the 
corresponding mean initial conditions a™^"" and a^^""^ can be significantly differ- 
ent, due to different initial mean longitudes (or similarly different initial resonant 
angle values). 

it can not induce a displacement of the equilibrium points in the phase space. 
Actually, a clever explanation can be found regarding the way the sampling 
is considered in the elaboration of the graphics. More specifically, it is worth 
noting that, at first, the sampling is carried out with respect to osculating 
initial conditions. 



Second, within the framework of mean motion theory, it is well-known that, 
due to the short-period oscillations, the mean and the osculating initial con- 
ditions can not be considered to be equal. In other words, for the same fixed 
value of the initial osculating semi-major axis and for various initial mean lon- 
gitude, we obtain different values for the mean semi-major axis; as explained 
with Figure [71 Actually, the different initial mean longitudes induce a phase 
difference in the corresponding evolution of the semi-major axis, leading to 
different mean initial semi-major axes. Let us remark that the maximum dif- 
ference between both the mean semi-major axes is directly related to the order 
of magnitude of the short-period variations, and, as a consequence, is also di- 
rectly related to the area-to-mas ratio. 



More rigorously, the difference between osculating and mean initial conditions 
is a well-defined transformation, depending on the generating function used 
within the averaging process allowing to change from mean to osculating dy- 
namics. For further details conc erning t his ex plicit transforination, we refer to 
the Lie algorithm discussed in iDepritI (119691 ) and iHenrardI (jl970l ). However, 
because we bound our analysis mainly to numerical simulations, we cannot 
access such generating function; we can nevertheless overcome this problem 
by numerically computing, for each semi-major axis osculating initial condi- 
tion, the related mean initial semi-major axis, by considering the average over 
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Fig. 8. Relation between the mean semi-major axis and the resonant angle for various 
values of the osculating semi-major axis. The first osculating semi-major axis is 
taken above the libration region, the second is related to an osculating semi-major 
axis sampling which crosses the libration region and, finally, the third sampling is 
taken below this region. 



a short time span of 10 days. As an illustration, in Figure [H we give the 
relation between the mean semi-major axis and the resonant angle for var- 
ious values of the osculating semi-major axis {A/m = 10 m^/kg). The first 
difference is related to a semi-major axis sampling taken above the libration 
region, the second is related to a semi-major axis sampling which crosses the 
libration region and finally, the third sampling is taken below this region. In 
conclusion, we clearly see that the order of magnitude of the differences is, 
as previously mentioned, the order of the amplitudes of the daily variations 
observed in the semi-major axis dynamics. Let us note that in the latter case, 
i.e. A/m = 10 m^/kg, the differences reach at most 27 km, which correspond 
exactly to the difference between the stable and unstable equilibrium points, 
as shown in Figure [6] (bottom, left). 

We can thus apply numerically the transformation as a post-treatment process, 
that is considering the MEGNO values not in the osculating initial conditions 
phase space, but in the mean initial conditions phase space. For the sake of 
comparison with Figure [6l we show the results once such a transformation has 
been applied (Figure [9]): it is clear that now the vertical gaps between both the 
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Fig. 9. The MEGNO computed as a function of initial mean longitudes Aq and initial 
mean semi-major axis qq. The model is the same as in Figure [6l The area-to-mass 
ratio is A/m = 5 and 10 m^/kg for the left and for the right graph, respectively. 

stable and unstable equilibrium points are almost completely eliminated, hence 
these points have almost the same mean semi-major axis, getting rid of the 
what we called the ^^short-period artefacf . The thin light waves crossing the 
Figure [9] are due to gaps in the set of initial conditions and have no dynamical 
significance (also valid for the only one light wave crossing the Figure [13] in 
section [6]). Let us also remark that, from now on, all the results will be shown 
in the mean initial conditions phase space. 



5. 3 Initial time at epoch and importance of the mean eccentricity 



One should also recall that solar radiation pressure leads to a theoretical 
equilibrium defined both in eccentricity cq and longitude of perigee wq. The 
conditions leading to such an equilibrium can be written as 

3 A 1 e A 

Co = -CrPr COS^ - ^ 0.01 Cr — , 

2 m nariQ 2 m 

where n and Uq are the angular motions of both the space debris and the Sun 
respectively, e is the obliquity of the Earth with respect to the ecliptic and 
A0(O) the initial ecliptic longitude of the Sun. If these conditions are fulfilled, 
it has been shown (Chao, 2006, and later Valk et al., 2008), that the eccen- 
tricity vector (e cosci7,e sinzu) remains constant, leading to a fixed value of 
both the eccentricity and longitude of perigee. 



As an illustration. Figure [TO] shows the mid-term variations of the eccentric- 
ity for a fixed value of the area-to-mas ratio {A/m = 10 m^/kg) and fixed 
initial conditions, namely, ao = 42164 km, cq = 0.1, io = rad, Qq = ujo = 
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Ao = rad. It is clear that, apart from a phase difference, the amphtudes 
of the variations of the eccentricities are quahtatively the same, except when 
adopting an initial time at epoch equal to 21 March. In this latter case, the 
eccentricity remains almost constant, as expected by the theory. 
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Fig. 10. Mid-term variations (yearly oscillations) of the eccentricity for fixed ini- 
tial conditions ao = 42164 km, eo = 0.1, = rad, ^Iq = loq = Xq = rad and 
A/m = 10 m^/kg. Various initial times at epoch to, namely different initial eclip- 
tic longitudes of the Sun A0(O), were used for the numerical propagations. Only 
J2, C22, 5*22 and direct solar radiation pressure taken into account. 



Figure [TT] shows the phase space in mean semi-major axis and longitude for 
A/m = 10 m^/kg and fixed values of the initial conditions, namely Cq = 0.1, 
io = 0.004 rad, Qq = Uq = rad. The differences between the two graphs only 
depends on the initial time at epoch parameter to- We could actually expect 
that different initial times at epoch, namely, different initial ecliptic longitudes 
of the Sun A0(O), will reveal a quite rich coUection of behaviors, depending on 
the different states with respect to the before-mentioned eccentricity equilib- 
rium. Actually, assuming an initial time at epoch of 21 December 2001, we see 
clearly that the phase space is filled by a large number of chaotic orbits (Fig- 
uredH left). On the contrary, starting with an initial time at epoch of 21 March 
2000, that is adopting a Sun pointing longitude of perigee (A0(O) = rad), 
the MEGNO values tend to be smaller and associated with significantly narrower 
chaotic regions, always located close to the separatrices (Figure [11], right). In 
the latter case, the eccentricity presents only small yearly variations due to 
the proximity of the theoretical equilibrium. 
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Fig. 11. MEGND computed as a function of initial mean longitude Aq and semi-major 
axis oq. The equations of motion include the central body attraction, the second 
degree and order harmonics J2 , C22 and S22 , the luni-solar interaction as well as the 
perturbing effects of solar radiation pressure. The mean longitude grid is 1° and 
the semi-major axis grid is 500 m spanning the 42164 it 35 km range. The initial 
conditions are cq = 0.1, io = 0.004 rad, = ujq = Q rad with A/m = 10 m^/kg. 
The patterns have been obtained using two different initial times at epoch, namely 
21 December 2000 (left) and 21 March 2000 (right), respectively. 

Therefore, these results seem to suggest that high amplitude variations of the 
eccentricity increase considerably the extension of chaotic regions close to the 
separatrices and, conversely, small eccentricity variations seem to minimize 
considerably the extent of chaotic regions. To justify this assumption, we per- 
formed a dedicated numerical simulation with the same set of parameters used 
in the one reported in Figure [TT|, but considering higher values of the initial 
eccentricity. The results are reported in Figure [121 the chosen time at epoch is 
21 December 2000 and the initial eccentricities are, eo = 0.2 (left panel) and 
Co = 0.4 (right panel). In the latter case, the huge variations of the perigee 
altitude, induced by the large variations of the eccentricity as well as by the 
variations of the semi-major axis, leads to even more complicated dynamics. 
These results thus confirm the importance of the initial eccentricity in the 
appearance of chaos. 



6 Secondary resonances 



It is worth noting that inspecting Figures [HI [ID and [T2] we clearly note the 
presence of some additional patterns located on both sides of the separatrices 
in the phase space. These never seen before regions, hence unexplained so far, 
are actually characterized by very low MEGNO values. Indeed, this observation 
underlines the fact that the dynamics of high area-to-mass ratio space debris 
is even more intricate than expected. In the following two sections we will 
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Fig. 12. MEGND computed as a function of initial mean longitude Aq and semi-major 
axis oq. The equations of motion include the central body attraction, the second 
degree and order harmonics J2 , C22 and S22 , the luni-solar interaction as well as the 
perturbing effects of solar radiation pressure. The mean longitude grid is 1° and 
the semi-major axis grid is 500 m spanning the 42164 it 35 km range. The initial 
conditions are io = 0.004 rad, Q.q = ujq = Q rad with A/m = 10 m^/kg. Time 
at epoch is 21 December 2000. The patterns have been obtained using two initial 
eccentricities, eo = 0.2 (left) and eo = 0.4 (right). 

provide some numerical results and an analytical theory, based on a simplified 
model, to better understand such zones. 



6.1 Numerical investigations 

We followed a large set of near-geosynchronous space debris, related to an 
extremely large set of initial conditions taken on both sides of the pendulum- 
like pattern, and for each one of the 72 000 orbits we computed the related 
MEGNO indicator. The initial conditions have been fixed by a mean longitude 
grid of 1°, spanning 360°, and a semi-major axis grid of 1 km, spanning the 
42 164±100 km range, while the remaining orbit parameters and time at epoch 
are the same as in Figure [61 Moreover, as in the previous extended analyses, 
the model of forces also includes the central body attraction, the second de- 
gree and order harmonics J2, C22 and S22 as well as the combined attractions 
of the Sun and the Moon. The perturbing effects of direct solar radiation pres- 
sure are also taken into account for a high area-to-mass ratio fixed at 10 m^/kg. 

The results are reported in Figure [T3l which is nothing but an extensive en- 
largement of the phase space presented in Figure [6] (bottom, left). This phase 
space widening clearly underlines the before-mentioned additional structures 
located at ± 40 km on each side of the resonant area. Furthermore, besides 
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these patterns, what is of special interest is that this Figure also brings to 
the light supplementary structures located at approximately 80 km on both 
sides of the main resonance, suggesting that the phase space is actually fo- 
liated by a larger set of secondary structures. Moreover, the width of these 
additional patterns and the numerical values of the MEGNO both seem to be 
directly related to the inverse of the distance with respect to the resonant area. 




Resonant angle [degree] 



Fig. 13. MEGND computed as a function of initial mean longitude Aq and semi-major 
axis qq. The equations of motion include the central body attraction, the second 
degree and order harmonics J2, C22 and S22 as well as the luni-solar perturbations. 
The mean longitude grid is 1° and the semi-major axis grid is 1 km, spanning 
the 42164 it 100 km range. The initial conditions are cq = 0.002, iq = 0.004 rad 
and r^o = '^o = rad. The area-to-mass ratio is 10 m^/kg. Time at epoch is 
25 January 1991. 



In addition, we also performed a set of similar numerical investigations, in or- 
der to distinguish qualitatively the relative relevance of some parameters such 
as the initial mean eccentricity, the value of the area-to-mass ratio, as well 
as the importance of the 1:1 resonance and of the third-body perturbations 
in the occurrence of such secondary structures. Even though these results are 
not presented here in detail, we can draw the following preliminary conclu- 
sions: the second order harmonic J2, as well as the third-body perturbations, 
do not seem to be really relevant and crucial in the appearance of these ad- 
ditional patterns. In other words, the unexpected patterns occur only when 
taking into account the combined effects of both the second order and degree 
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Fig. 14. Blow-up of the phase space with the specification of a resonant angle sec- 
tion (horizontal black solid line), that is the set of orbits having the same (oscu- 
lating) initial resonant angle value, near the first stable equilibrium point, namely 
^section _ gi Qjo ^^^^ panel). Evolution of MEGNO with respect to the initial semi- 
major axis oo for the specified section (middle panel). The fundamental period of 
cr with respect to the initial semi-major axis ao, computed by means of frequency 
analysis for the specified section (bottom panel). The estimation of the periods are 
made over a 20 years period of time. 

harmonic and direct radiation pressure. As a matter of fact, the extended 
numerica l investigations perf ormed in Figure [6] (top, left), or similarly those 
shown in Breiter et al. (l2005h . also present these structures, even though they 



are difficult to perceive. Actually, the extension and chaoticity indicator of the 
secondary patterns seem to be directly proportional to the area-to-mass ratio 
value or, equivalently directly proportional to the mean value of the eccentric- 
ity. 



To get even more concluding results, we considered a blow-up of the phase 
space (dashed line rectangle in Figure [T5]) with really high resolution sam- 
pling (150 meters in the semi-major axis a and 0.3° in the resonant angle a). 
Figure [14] (top) shows this phase space widening wherein we defined a so- 
called resonant angle section (horizontal black solid line), that is the subset 
of orbits having the same initial resonant angle value. This resonant angle 
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section spans the complete range in semi-major axis and passes close to the 
stable equilibrium point. For each orbit defined on this section, we computed 
the MEGNO indicator and in Figure [H] (middle) we report this value at the end 
of the simulation as a function of the semi-major axis. 

To double check our results, we performed a frequency analysis investigation 
(see Laskar, 1990 and 1995, and Noyelles et al., 2008) aimed to study the 
behavior of the proper frequency of the resonant angle a, whose results are 
reported in Figure [H] (bottom). Here one can clearly notice the distinctive 
characteristics regarding the well-known 1:1 resonance between the mean lon- 
gitude and the sidereal time. Indeed, both MEGNO and the fundamental period 
show distinctively a minimum close to the stable equilibrium point. In this 
case, as previously mentioned in Section HJ MEGNO should slowly converge to 
Y{t) = 2 everywhere, except at the equilibrium point where the limit value is 
Y{t) = 0, that's why, using a finite integration time, we obtain such V-shaped 
curve, close to in the center of the resonance and to 2 on the borders. It is 
also worth noting that the fundamental period of a is reported to be close to 
2.25 years, which is in good agreement with the well-known 818 days libration 
period of a typical uncontrolled near-geosynchronous object. Near the sepa- 
ratrices, MEGNO clearly presents some obvious high values which confirms the 
presence of chaotic orbits. Here, the fundamental period reaches significant 
values and, as a matter of fact, is not well determined, once again supporting 
the result of the existence of a chaotic zone. 



Moreover, the use of frequency analysis allows us to support strongly the 
hypothesis that the additional patterns are actually related to secondary 
resonances . Indeed, if we look at the evolution of the fundamental period 
with respect to the semi-major axis, it is clear that the so-called secondary 
resonances are associated, regarding the angle a, with periods which are com- 
mensurate with 1 year. More precisely, the major secondary resonances, lo- 
cated at approximately 40 km on both sides of the pendulum-like pattern, are 
related to a 2 years fundamental period of a. Concerning the farther patterns 
located at ±80 km, the fundamental period of a turns out to be very close to 
1 year. As a consequence, we can presumably assume that these secondary res- 
onances are actually related to a commensurability between a and the 1 year 
period angle Xq, that is the ecliptic longitude of the Sun. 

To justify this assumption, we focused our attention to the major secondary 
resonances located at ±40 km on both sides of the pendulum-like pattern, 
considering the time evolution of various linear combinations of a and A0. 
For this purpose, we considered various initial semi-major axes in the phase 
space. The results are shown in Figure [151 At first glance, it is apparent that 
three propagations stand apart from the others. In the first row of Figure [151 
that is regarding the evolution of the resonant angle a, we clearly identify the 



28 




10 20 30 10 20 30 10 20 30 10 20 30 10 20 30 



Time [years] 

Fig. 15. Time evolution of the angles a, 2cj + A0 and 2a — Xq (in radians) for several 
semi-major axes. In the lower major secondary resonance, oq = 42 122 km. In the 
eye of the principal resonance, ao = 42 188 km. Between the primary resonance 
and the upper secondary resonance, oq = 42 204 km. Inside the upper major sec- 
ondary resonance, oq = 42 212 km. Outside the upper major secondary resonance, 
ao = 42 230 km. 

well-known characteristics related to the primary resonance. In particular, in 
Figure [T5k. that is when considering an initial semi-major axis inside the pri- 
mary resonant (oq =42 188 km), a shows the well-known long-period libration 
(2.25 years), whereas a circulates outside this region. Furthermore, what is 
of special interest is the time evolution of both 2cr + A0 and 2cr — A0, shown 
in the second and third row, respectively. It is clear that most of the time 
these angles show a circulation regime. However, when considering an initial 
semi-major axis inside the major lower secondary resonance for 2a — Xq or, 
similarly inside the major upper secondary resonance for 2a + Xq, both these 
angles show a significant long-term evolution (Figure [T5b and [T5b) . 



6.2 Analytical investigation - simplified model 



The presence and the location of these secondary resonances can be studied 
using an appropriate simplified model. Hence we model the averaged geosta- 
tiona ry motion by a pen dulum-like system, given by its Hamiltonian formula- 
tion (jValk et al.l . |2008q) up to order in the series expansion 



n 



2L2 



-0L 
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where 



fia and S'22oo(^; ^) = ^^22 cos 2a + S'22 sin 2cr . 



In the context of direct solar radiation pressure, we can introduce the factor 
Z proportional to A/m through the eccentricity e (for further details, we refer 
to the averaged simplified analytical model developed in Valk et al., 2008). As 
a first approximation, the time evolution of both the eccentricity e and the 
longitude of perigee zu were found to be (neglecting the obliquity of the Earth 
with respect to the ecliptic) 

Z 

e cos zu = cos Aq + ao , 

LriQ 

Z 

e sin w = sin Xq — (3o , 

LUq 

which introduces A0 in the Hamiltonian. The quantity Hq is the mean motion 
of the Sun and both ao and Pq are related to the initial conditions with respect 
to the eccentricity and the longitude of perigee). The resulting Hamiltonian 
takes the generic form 

n = -eL + — cos (2a - 2ao) - — 2 cos (2a - 2ao) cos (Aq + 6) , 

where 6, F, G, ao are constants. A suitable transformation is then necessary to 
introduce action-angle variables {ip, J) in the libration and in the circulation 
region of the double pendulum, in such a way any trajectory of the double 
pendulum is characterized by a constant action J and a corresponding constant 
frequency ip- Rewriting the perturbed system (because of A0 terms) by means 
of these new variables and then using the expansions in Bessel functions, we 
could isolate any resonance of the type kip ± A© in the circulation region, 
for any \k\, and in the libration region, for \k\ > 3, which corresponds to 
our frequency analysis. This analysis is surely promising, but it is outside the 
goals of this paper. Further investigations will be det ailed in a forthcoming 
publication ( Lemaitre et al. . submitted for publication! ). 



7 Conclusions 



The predictability of the trajectory high area-to-mass ratio space debris lo- 
cated near the geosynchronous region was investigated by means of a recent 
variational chaos indicator called MEGNO. Thanks to this technique, we clearly 
identified the regular (stable) and irregular (chaotic) orbits. This efficient 
method allowed us to obtain a clear picture of the phase space, hence show- 
ing that chaotic regions can be particularly relevant, especially for very high 
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area-to-mass ratio objects. Moreover, we discussed the importance of both the 
initial eccentricity and time at epoch in the appearance of chaos. 

Finally, we brought to the fore a relevant class of additional unexpected pat- 
terns which were identified as secondary resonances, that were numerically 
studied by means of both the MEGNO criterion and frequency map analysis, 
to eventually conclude that they involve commensurabilities between the pri- 
mary resonant angle and the ecliptic longitude of the Sun. We also presented 
an analytical scheme that could explain their existence. It will be the subject 
of further work. 
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